Stress distribution and the fragility of supercooled melts 
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We formulate a minimal ansatz for local stress distribution in a solid that includes the possibility 
of strongly anharmonic short-length motions. We discover a broken-symmetry metastable phase 
that exhibits an aperiodic, frozen-in stress distribution. This aperiodic metastable phase is char- 
acterized by many distinct, nearly degenerate configurations. The activated transitions between 
the configurations are mapped onto the dynamics of a long range classical Heisenberg model with 
6-component spins and anisotropic couplings. We argue the metastable phase corresponds to a 
deeply supercooled non-polymeric, non-metallic liquid, and further establish an order parameter 
for the glass-to-crystal transition. The spin model itself exhibits a continuous range of behaviors 
between two limits corresponding to frozen-in shear and uniform compression/dilation respectively. 
The two regimes are separated by a continuous transition controlled by the anisotropy in the spin- 
spin interaction, which is directly related to the Poisson ratio a of the material. The latter ratio 
and the ultra-violet cutoff of the theory determine the liquid configurational entropy. Our results 
suggest that liquid's fragility depends on the Poisson ratio in a non-monotonic way. The present 
ansatz provides a microscopic framework for computing the configurational entropy and relaxational 
spectrum of specific substances. 



If cooled sufficiently rapidly, a liquid may fail to crystal- 
lize, but will instead remain in a metastable, supercooled 
state. Upon further cooling, the relaxation times in a su- 
percooled liquid grow very rapidly as the mass transport 
becomes activated, in contrast with the mainly collisional 
transport near the fusion temperature. Because the lo- 
cal structures are much longer-lived than the vibrational 
equilibration times, the activated-transport regime repre- 
sents a state with a broken translational symmetry, even 
though the corresponding, aperiodic structure shows no 
obvious distinction from a snapshot of an ordinary, uni- 
form liquid. (It is said the heterogeneity is "dynamical.") 
Since the symmetry is broken gradually with lowering the 
temperature - beginning with the highest frequency mo- 
tions - a transition into this "aperiodic-crystal" state is 
not sharp, but, instead, is a soft cross-over centered at 
a temperature T CI ^^ corresponding universally to vis- 
cosity 10 Ps or so.-i 3 - The cross-over into the activated 
regime is a finite-dimensional analog of a mean-field ki- 
netic catastrophe of the mode-coupling theory (MCT), 
whereby the motional barriers would diverge at a tem- 
perature Ta, even though the configurational entropy is 
still extensivei^ 

In additional contrast to the mean-field transition at 
Ta, the cross-over at T cr exhibits two emerging length 
scales: One length scale is the molecular length a that 
signifies the volumetric size of a chemically rigid unit - of- 
ten called the "bead" - that is not significantly perturbed 
during activated transport. Conversely, the beads inter- 
act with each other weakly, comparably to the Lennard- 
Jones interaction.— The bead may be thought of as a 
coarse-graining length, beyond which activated motions 
are largely independent of chemical detail, but fully char- 
acterized by a single, bulk quantity. This bulk quantity 
is the excess liquid entropy relative to the corresponding 
crystal, usually called the configurational entropy. The 
magnitude of the configurational entropy per bead, s c , 
directly gives the number of alternative aperiodic con- 



figurations available to a region of a supercooled liquid 
containing N beads, i.e. e ScN ^ ks . The bead usually con- 
tains two-three atoms, but could be bigger for molecular 
liquids containing large rigid units such as benzene. The 
other length scale emerging during the crossover is the 
so called Lindemann length d^,— which is the molec- 
ular displacement at the mechanical stability edge. This 
length is nearly universal: dh ~ a/10, and characterizes 
bead displacements during transitions between distinct 
aperiodic packings in the metastable, aperiodic crystal 
phase. One may view the crossover into the activated 
liquid regime as a "localization" transition, whereby the 
emerging metastability of local structures is signaled by a 
discontinuous transition from a uniform liquid to a state 
with a non-zero force constant of the Einstein-oscillator. 1 
The random first order transition (RFOT) theory utilizes 
this view to analyze the activated transport in chemically 
distinct (non-polymeric) fluids in a unified fashion^— 
(See Ref.— for a review.) The RFOT theory predicts that 
in a fully developed activated regime, the structural re- 
laxation time is determined solely by the configurational 
entropy per beadi^ 
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where T v ib r — 1 psec is the microscopic time scale char- 
acterizing vibrational relaxation. By Eq.([T]), system- 
specific deviations from the pure Arrhcnius temperature 
dependence of r result from variations in the value of 
the heat capacity jump at the glass transition tempera- 
ture T g per bead: Ac p = T(ds c /dT)\ T = Tg - The so called 
fragility index m = d(logr)/Td(l/T)|x=T g gives a quan- 
titative measure of that deviation. Small and large de- 
viations from the Arrhenius T-dependence (correspond- 
ing to small and large m respectively) are often called 
strong and fragile behaviors . 12 ' 13 The RFOT theory pre- 
dicts m ~ ;i 1.7Ar ; ,.- fl 1 1 in excellent agreement with ex- 
periment. 

Hall and Wolynes have put forth a simple specific 



2 



model that relates the degree of molecular connectivity to 
the heat capacity jump Ac p . 15 Establishing such connec- 
tions between local-chemistry and thermodynamics for 
actual substances is difficult, however, hampering our ef- 
forts to make first principles estimates of the configura- 
tional entropy and, ultimately, the glass-forming ability 
of those substances. Indeed, since the viscosity is directly 
related to the average relaxation time,— 
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the crystal nucleus growth, which is subject to viscous 
drag^i is determined by the configurational entropy s c 
per bead. The actual chemical identification of beads 
is relatively straightforward in molecular, but less so in 
covalently bonded compounds, where, for instance, the 
apparent bead often corresponds to a non-integer frac- 
tion of a stoichiometric unitj^^ We note that these chal- 
lenges are, of course, not unique to theories of the glass 
transition but equally pertain to quantitative descriptions 
of the ordinary liquid-to-crystal transition. Yet perhaps 
in reflection of these open questions, which we view as 
quantitative, many believe fragile and strong behaviors 
actually have distinct mechanisms. 

In an effort to develop a microscopic description of the 
configurational dynamics in specific substances, here we 
propose and work out several consequences of a novel 
quasi-continuum ansatz for local stress distribution that 
implements direct interactions in the activated regime 
semi-empirically, via the local elastic properties of the 
material. The ansatz incorporates the possibility of local, 
short-length motions that are similar in spirit to Einstein 
oscillators, but are strongly anharmonic and account ex- 
plicitly for the tensorial nature of the relation between 
local deformation and stress in solids. In particular, ex- 
plicitly treated is (high- frequency) shear resistance, which 
is characteristic of deeply supercooled liquids and, gener- 
ally, activated liquid transport. This fully tensorial treat- 
ment of stress may be viewed, among other things, as a 
systematic improvement on self-consistent scalar-phonon 
theories of aperiodic crystalsj^^ 

The many-body effects that lead to the emergence of 
shear resistance in liquids at high densities do not lend 
themselves easily to perturbative treatments that use the 
uniform liquid as the reference state. In contrast, we use 
a fully mechanically stable state as the reference state, 
and then uniformly allow for local, short-length anhar- 
monic motions. Within this ansatz, a metastable, ape- 
riodic frozen-in stress pattern emerges self-consistcntly, 
whereby the bead size is identified as the typical length 
scale of the stress heterogeneity. This result implies 
that even in a covalently bonded material, the concept 
of a bead is entirely unambiguous. The frozen-in stress 
pattern is multiply degenerate, whereby the transitions 
between alternative configurations can be mapped onto 
the dynamics of a classical Heisenberg model with six- 
component spins. The six components correspond to 
the independent entries in the local frozen-in deforma- 
tion tensor. The average length of the spin emerges as 



an order parameter for the aperiodic-to-periodic crystal 
transition. The transitions between different spin config- 
urations are rare events implying the corresponding liq- 
uid states are long-lived and thus supporting the view 
of liquids in the activated-transport regime as aperiodic 
crystals. In other words, we explicitly confirm the non- 
trivial notion of the RFOT theory that aperiodic molec- 
ular assemblies can exhibit multiple states and, at the 
same time, support shear. 

Already in a mean-field analysis of the spin model, we 
establish that the broken-symmetry state can exhibit a 
range of behaviors interpolating between frozen-in shear 
and uniform compression/dilation. These two limits cor- 
respond to a 5-component Heisenberg- and Ising-like fer- 
romagnets respectively. The resulting heat capacity vari- 
ations, together with the RFOT-derived intrinsic relation 
between thermodynamics and kinetics, as reflected in the 
to vs. Ac p relation, implies that the spin model exhibits 
a broad range of fragile-to-strong behaviors during acti- 
vated rearrangements. The heat-capacity jump per bead 
of the corresponding liquid will also depend on the ul- 
traviolet cutoff of the theory, as specified by the bead 
size. We will observe that the fragility depends on the 
Poisson ratio of the material, but in a complicated, non- 
monotonic fashion. In any case, in view of the appar- 
ent broad range of strong-to-fragile behaviors exhibited 
in the present model, the present analysis strongly sug- 
gests the activated dynamics have the same mechanism in 
both strong and fragile liquids, consistent with the basic 
picture of the RFOT theory. 

A potential formal benefit of the present approach is 
that it is possible to map the activated-transport liquid 
regime onto a spin model on a fixed lattice, thus opening 
the possibility of applying to liquids the many computa- 
tional techniques developed in the context of spin models. 
Importantly, this mapping may provide a first principles 
basis for computing the configurational entropy for spe- 
cific substances and the rate of crystallization from the 
supercooled state, since elastic constants at the length 
scale of a bead, which are the input into the theory, can 
be determined by ab initio methods. 



I. MICROSCOPIC ANSATZ FOR LOCAL 
STRESS DISTRIBUTION 

Informal motivation. Consider a substance that can 
both crystallize and vitrify and assume for concreteness 
that the glass has a greater specific volume than the corre- 
sponding crystal, as is usually the case. A quenched melt 
or frozen glass of such a substance may be thus viewed as 
a result of the following action: Apply negative pressure 
to slowly expand the crystal until the lattice relaxes into 
one of the myriad available glass structures, and then re- 
lease the negative pressure. (We assume cracking does 
not occur.) The sample will remain in the metastable, 
glass state; the eventual transition to the lower free en- 
ergy, crystal state will occur by nucleation and is sub- 
ject to the surface tension between the two phases. In 
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FIG. 1: Illustration of the notion of an aperiodic metastable 
state as a dilated crystal that was allowed to relax. 



a standard fashion, the surface tension implies that the 
dependence of the free energy on the extent of dilation 
exhibits a non-concave portion^ (By "concave," we will 
consistently mean concave- up.) 

We can see qualitatively how the non-concavity arises 
by estimating the free energy of the solid using a 
simple non-linear Einstein-oscillator model. Suppose 
the configurational partition function of each oscilla- 
tor is d 3 re~ kr / 2kBT , where we incorporate the non- 
linearity qualitatively by assuming a hard lower cut-off 
in the configurational integral, at a displacement g. k 
is the force constant in the harmonic limit. Denote the 
force constants in the crystal and the expanded state as 
fcx and fee respectively. Clearly, Uq < fcx. We sketch 
the free energies of the states as functions of g in Fig|TJ 
taking into account that the curve corresponding to the 
expanded state is shifted upward from the crystal state 
by the free energy cost of expansion. According to Fig|T] 
under sufficient negative pressure (essentially equal to the 
negative slope of the common tangent of the two curves), 
the crystal will indeed begin converting into a collection of 
non-linear oscillators. The vast majority of the configura- 
tions of these oscillators are of course aperiodic, implying 
the metastable phase is aperiodic too. This qualitative 
free energy graph tentatively confirms the view above of 
a quenched-melt/frozen- glass as a sufficiently expanded 
crystal that was allowed to relax. The graph also sug- 
gests an appropriate order parameter for the transition, 
namely the amplitude of short-wavelength, strongly an- 
harmonic motions. 

In writing a formal ansatz for a solid that allows 
short length anharmonic motions, we first point out 
that local mechanical instabilities are perfectly compat- 
ible even with macroscopic rigidity, examples includ- 
ing orientational glasses^ and ferroelectrics of order- 
disorder typej22~— The latter systems consist of rigid, 
periodic scaffolds where a subset of local motions could 
be marginally stable (resulting in a vanishing frequency 
of transverse optical phonons ? -) or even multistable. In 
contrast, we wish to consider here a situation where all 



short length motions are allowed to be multistable. The 
resulting material may or may not be macroscopically 
stable, as will be determined self-consistently. We will 
limit ourselves to a dynamical regime, in which relative 
distortions of individual bonds during lattice reconfigura- 
tions are only within 10% or so^i as appropriate both 
for activated transport and stable crystals. Such small 
distortions imply the harmonic approximation for local 
stress is semi-quantitative. Take the Landau local free en- 
ergy density as a function of the local deformation profile 
Uif f = \uij A-ijki Uki, where Ayiy is the standard elas- 
tic tensor and the Einstein's summation convention for 
Latin indeces is implied^ (utj = \{pUij dxj + duj/dxi), 
i, j = 1, 2, 3, where u is the local displacement relative to 
a steady-state reference state.) 

Decompose the deformation tensor into long- and 
short- wave length components: u%j = ufj + <iy, where 
uf- only includes Fourier components with wave-vectors 
within a certain cut-off surface in fc-space, while dij con- 
tains only Fourier components outside of the cut-off sur- 
face. Let us evaluate the free energy of the system sub- 
ject to a constraint on the self-energy of the local motions 
dij A-ijkl dki = 9 , thus forcing the short-length motions to 
be essentially anharmonic Einstein oscillators. If the cut- 
off surface in the Fourier space is generically at k ~ fc m icro, 
the anharmonic oscillators fill the space, approximately 
one per volume (7r/fc m j cro ) 3 , but not necessarily in a peri- 
odic fashion. The detailed lattice in which the oscillators 
are arranged depends on the precise shape of the cut-off 
surface. Since the presence of strongly anharmonic de- 
grees of freedom generally requires expanding the lattice 
from the reference state, owing to stcric repulsion, we 
must include in the full free energy of our solid a penalty 
for increasing g, which we denote as F ex (g). The pre- 
cise form of this function is not essential for the majority 
of our conclusions; clearly it is a monotonically increas- 
ing, concave function of g, to ensure the compressibility 
is positive. As a result, the full free energy is: 

F(g) = F (g) + F^g), (3) 

where the free energy of the ansatz Fq is computed by 
summing over all configurations of local stress at temper- 
ature T = l/k B (3: 

e -PF (g) = e -0f*Vf(T)^ ( 4 ) 

{u(r),d(r)} 

subject to the constraint 

dij(r)Aij k idki{r) = g 2 {r). (5) 

The local free energy density f(r) of a particular stress 
configuration at site r is 

f( r ) = 7^( u ij +dij)Aijki {u k i +d k i). (6) 

We have dropped the superscript "<" from the long- 
wavelength, "acoustic" component of the total deforma- 
tion u < , to simplify notation. For the sake of argument, 
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we will assume the order parameter is uniform: 

9\r)=g\ (7) 

Formal extension to a non-uniform g[r) is straightfor- 
ward. 

In the following, it will be often convenient to present 
the six independent entries of the symmetric tensors u,j 
and dij as 6-component vectors, which is sometimes re- 
ferred to as the Voigt (or Voigt-Mandel) notation ! 26 ' 27 
Here we choose the specific realization of the Voigt nota- 
tion, in which the tensor A is a second rank (positively de- 
fined) tensor of size 6x6. Upon a suitable linear transfor- 
mation, see Appendix [3] we can formally rewrite Eq.© 
as 



f(r) 



■9) 



where 



E 



</> 2 = Uij Aijkl Ukl, 



(8) 



(9) 



and 



Hal 



= dij Aijki dki- 



(10) 



Throughout the article, we will consistently indicate 6- 
component vectors with arrows: g = Y] a — i 9ae a , where 
e a S/3 = 5 a p; and three component vectors with bold sym- 
bols: u — Uiei, where e^ej = 5ij. We will not apply the 
Einstein summation convention to Greek indeces. 

Even though we have written g(r) (or (kj{r)) as a con- 
tinuum field, it is understood that g(r) actually stands 
for a discrete vector located at the site of the lattice of 
the non-linear Einstein oscillators that is closest to the 
point r in space. Switching between discrete summation 
over lattice sites and continuous integration over space is 
straightforward: If the volume occupied by a single vector 
g is a 3 - implying fc m i cro — 7r/a - we may interchange the 
summation and integration according to a 3 J] f> J dV. 
Since vectors g correspond to 6 degrees of freedom, the re- 
gion corresponding to one vector g must contain at least 
two atoms. (f>(r) (or iiy(r)) is of course, too, an ap- 
proximation to acoustic modes of a discrete lattice, at 
wave-lengths exceeding 2ir/k m i cro . Note that the number 
of these acoustic modes is three per vector g(r), corre- 
sponding to the three translational degrees of freedom of 
the vector-containing region. 

For the sake of argument we will consider the simplest 
case of an isotropic reference state: 



Aijki = \6ij$kl + fJ>(6ik$ji + SuSjk), 



(11) 



where A and fj, are the Lame coefficients.— /i is also called 
the shear modulus. In the isotropic case, the components 
of vectors <f> and g have a particularly lucid meaning. 



One component - called "dilatational" - corresponds to 
uniform compression: 

g\ = Kdl, (12) 

and the rest five - called "isochoric" - to pure shear: 

6 



(13) 



Q = 2 



where K = [A + (2/3)//] is the usual bulk modulus^ 



II. GLASS AS A FROZEN-IN, MULTIPLY 
DEGENERATE STRESS PATTERN 



We inquire about the possibility of (meta-stable) states 
with non-zero g by testing for non-concavity of F(g) from 
Eq. (J3]) . The latter peculiarity is only possible if F (g) has 
an inflection point. An analytic evaluation of -Fn(ff) ap- 
pears difficult (except in the mean- field limit, see next 
Section). Despite the simple bilinear form in Eq.®, the 
functional integration over u(r) in Eq.Q at a fixed con- 
figuration d^ (r) could not be performed by simply chang- 
ing variables u = u' — d and performing a Gaussian inte- 
gration with respect to the new variable u': In contrast to 
the purely harmonic, long- wavelength motions it, the new 
displacement field u' would now contain short-wavelength 
motions, which, by construction, are anharmonic. 

It will suffice for our purposes here to extract the small 
and large g asymptotics of F (g), which one may accom- 
plish by first averaging the Boltzmann weight f dV f( r ^ 
from Eq. @ over the directions of the g vectors on each 
site. By the aforementioned prescription for interchange 
between discrete summation over the oscillator sites and 
integration in space, we break up the integral in the 
exponent into a sum J2i a3 f( r i)- Note the formula 
(exp(xy)) = 7 m/2 _ 1 ( a ; 2 /)(xy/2) 1 - m / 2 r(m/2), where the 
averaging (...) is with respect to the mutual angle be- 
tween two m-component vectors x and y (see Ref.— or 
formula 9.6.18 of Refj^). Here, I„(x) is the modified 
Bessel function of the first kind of order v and T is 
the gamma function^ Using this formula, we average 
e~P a ^ over the directions of g on each site and obtain 
the following free energy at site r, per non-linear oscilla- 
tor: 



s 9 W = i[0 2 (r)+. 9 2 



J 2 [/3a 3 0(r). 9 ] 
Pa 3 "* [/3a 3 (j)(r)g] 2 ' 



In 



(14) 



up to an additive constant. The subscripts 0, g at 
f(p,g(f) above signify that the latter function is derived 
from the original free energy density f(r) from Eq.®, 
but now depends only on the absolute values of the 6- 
componcnt vectors <f> and g characterizing the long- and 
short-wavelength distortions respectively. We note the 
expression (1141) represents the exact free energy of our 
ansatz, up to the small numerical ambiguities related to 
the detailed shape of the cut-off surface in the Fourier 
space and the lattice of the spins. 
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To extract the small g behavior of Fo(g), we note that 
at small values of either field, expression (1T41) yields 



Ur) ^ l<f>\r) 
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(15) 



implying that the sole effect of the presence of the (small) 
short-wavelength anharmonic displacement of magnitude 
g is a renormalization of all acoustic frequencies down- 
wards by a factor of (1 - ff 2 /3a 3 /6) 1/2 ~ (1 - g 2 /3a 3 /12). 
Since the partition function of a classical harmonic os- 
cillator with frequency u) is (fcflT/fiw), this reduction 
in the acoustic frequency results in a free energy shift 
of — k B T(g 2 (3 a 3 / 12) per phonon. As mentioned earlier, 
the number of the acoustic modes is three times the 
number of vectors g, yielding that at small g, f(g) ~ 
g 2 /2 - 3{k B T/a 3 )(g 2 (3a 3 /12) = g 2 /4. Thus at small g, 
the second derivative of i*b(s) is positive. 

According to the asymptotic expansion of the Bessel 
function, I v (x) oc e x /x 1 / 2 ^ the two leading terms in the 
large g limit of Eq. (|14[) are given by 



UA r ) 



5k B T 



In g. 



(16) 



In contrast with Eq.®, the ^-dependent portion of the 
free energy is now uniform in space and thus no longer 
has short-wavelength components. As a result, we can 
always shift the u variable to eliminate the g depen- 
dence from the integrand in the partition function cor- 
responding to the free energy (fTo]) so as to make the ex- 
pression a pure quadratic in the new variable (f/: u = 
u + u' =>■ <f> = 00 + 0' where \<j> \ — g. (This lat- 
ter reference state physically corresponds to a uniform 
distribution of stress (not displacement!) at free energy 
density g 2 /2.) The resulting integration is <f>' becomes 
asymptotically Gaussian in the large g limit, despite the 
non-analytical dependence on u in Eq. (|16[) . As a result, 
fg(d) g ^°°> {5/2)(k B T/a 3 ) \ng, i.e. the second deriva- 
tive of Fo(g) becomes negative at large g, implying F${g) 
has an inflection point at an intermediate value of g. 

We thus conclude from the analysis of the asymptotic 
behavior of the free energy Fo(g) that with a suitable 
form of Fex(g), such as F cx (g) cx g 2 , it is possible for the 
full free energy F(g) to have a non-concave portion, also 
consistent with the mean-field solution in Section [ill! see 
Figf5] We reiterate that the presence of the non-concave 
portion means the system exhibits a metastable state at 
a non-zero value of g. The latter quantity therefore can 
used as an order parameter for a transition between the 
lowest free energy, fully stable state and the metastable 
state at non-zero g. 

The free energy of our ansatz, Eq. (fT4|) . helps under- 
stand the origin of the peculiar behavior of Fo(g). Indeed, 
at large enough magnitude of anharmonic displacement, 
g, the equilibrium value of the acoustic displacements be- 
comes non-zero. According to Eq. (Tl5)) this transition oc- 
curs at a critical value: 
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FIG. 2: Free energy profiles for local elastic strain at three 
representative values of the anharmonic displacement g, from 
Eq. (|14p . Above the critical value <? C r from Eq, (ll7|l . the lattice 
freezes in an aperiodic pattern at wave-length 2tt /fe m ; cro . 



We show the <fi dependence of the free energy from Eq. (114)1 
in Figl2l for several representative values of g. Note that 
the critical value of the excess free energy at a single 
vector, a 3 g 2 I /2 = 3k B T, is the total energy of a single 
atom in a harmonic lattice. 

In the metastable phase at g > 0, the thermodynamic 
potential from Eq.(|14j) is dominated by aperiodic config- 
urations of g 's, implying the modulation of the harmonic 
held is aperiodic, too. This potential becomes particu- 
larly vivid in the g 2 /k B T —> oo limit, i.e. (g — \<f>\) 2 /2 = 
[g — (uij Aijki Uki) 1 / 2 ] 2 /2. This form means that the ape- 
riodic phase corresponds to a frozen-in stress pattern at 
half- wave length a = 7r/fc m i cro , which we must thus asso- 
ciate with the bead size of the RFOT theory. Consistent 
with this identification, the form [g— (v,ij Ayti Uki) 1 ^ 2 ] 2 is 
similar to broken-symmetry functionals that arise in sys- 
tems with self-generated disorder ! 30 ' 31 Importantly, the 
free energy (| X4[) is an equilibrium quantity, implying the 
metastable aperiodic state corresponds to a melt above 
the glass transition. 

The aperiodic phase exhibits many distinct states sep- 
arated by activation barriers. To see this structural de- 
generacy explicitly, we first integrate out the acoustic de- 
formation in Eq.(j4|) at a fixed dij(r) configuration. Anal- 
ogous calculations have been performed, for instance, by 
Grannan at eh 32 ' 33 , who have considered defect-phonon 
couplings in the form of QijUij. Using their result or 
by straightforward path-integration in fc-space (see Ap- 
pendix [B}, we obtain the following effective Hamiltonian 
that couples the non-linear oscillators, in the tensorial 
and Voigt notations respectively: 



^2 a 3 K ijk i(r m - r„) d^d 



kl 



E E « 3 ^(*-m-r n )ffi OT) 4 n) ,(18) 

m<n a, /3=1 



(17) 

where the double (m < n) sums are over all bead pairs 
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and the spin-spin coupling is the Fourier transform: 

/" fc micro d 3 (ka) 

K ijk i(r)= / —i>coB{kr)K ijk i, (19) 

of the following tensor: 



flexible functional form for the non-linear oscillators, in- 
cluding higher order multipole terms in the spin-spin in- 
teraction etc., are possible, however will result in largely 
quantitative corrections. 



x + n 



„ _ I A A " 1 h~ A A 

^ijkl I I ^ijmp 1 ^klmq . l*-iim.n + *- 



\ + 2[l 



ijmp^klnq km 



(20) 

where k = k/k and the tensor h^u is from Eq. flTTT) . 
The Fourier image J a p(k) of the interaction J a p{r) from 
Eq. (|18[) is a dimensionless 6x6 matrix, which means here 
that the interaction between the spins scales as 1/r 3 at 
large distance, similarly to the usual dipole-dipole inter- 
action. 

We note that Hamiltonians similar to that in Eq. (fT5|) 
have been employed in the context of configurational 
glasses . 32 ! 34 Those treatments were restricted to defects 
of the type Qij cx [riirij — (l/3)<5y], where n is a three- 
component unit vector corresponding to the orientation 
of a polar group. Grannan at el* 3 - 2 - point out that mag- 
nets of the type in Eq. ((T8|) are frustrated even on periodic 
lattices. 

The classical Heisenberg-like ferromagnet from Eq. (fT5|) 
exhibits many distinct states separated by barriers, below 
its Curie temperature. (See next Section for the mean- 
field phase-diagram of the model.) These distinct states 
correspond to distinct configurations in the correspond- 
ing quenched liquid. The bead motions during transitions 
between these distinct states correspond to rotations of 
the six-component, classical "spins" g( m \ The relative 
displacement of two beads, say m and n, during a tran- 
sition is related to the transition-induced change in the 
vector (g ( m J — g^). The relation, however, is not obvi- 
ous because of the tensorial nature of the frozen-in stress 
resulting in vectors being six-component. As in any mag- 
net below its ordering transition, typical transitions be- 
tween distinct metastable states of the spin system from 
Eq. (TT5|) will involve the more spins the lower the tem- 
perature. Such bigger cooperative regions are necessary 
because the density of states for a region of fixed size 
decreases dramatically for lower energies. Specifically, 
according to the RFOT theory, the cooperativity size in 
the present case should increase as where s c is 

the entropy per bead (spin). Consistent with the coop- 
erativity, the individual vectors rotate little relatively to 
each other so as to minimize local strain, implying the 
beads will move nearly harmonically most of the time. 
The transition is, nevertheless, strongly anharmonic, as 
reflected in a high free energy barrier for the reconfigura- 
tion. An appropriate progress coordinate is the number 
of beads that have already moved ^ 

We conclude this Section by pointing out that the ob- 
tained mapping between activated liquid transport and 
the dynamics of the 6-spin Hamiltonian from Eq. (|18l) is 
not exact, of course, but nonetheless, does capture the es- 
sential features of the activated-transport regime in liq- 
uids. Systematic improvements, such as using a more 



III. MEAN-FIELD ANALYSIS OF THE SPIN 
MODEL 



It is straightforward to show that the Fourier image 
J a /3(k) of the interaction J Q ^(r) from Eq.(fT8]) has the 
following simple property: J 2 = J (i.e. p J aj3 J p-y — 

J Q7 ), and that the trace of the J matrix equals three for 
any value of the wave- vector k. These two notions imply 
the matrix J has two eigenvalues: 1 and 0, both triply 
degenerate. For comparison, the analogous matrix for 
the regular elastic dipole-dipole interaction between 3- 
component spins, kikj (see Appendix |B|) . has eigenvalues 
1, 0, 0. For the electric dipole-dipole interaction, the 
eigenvalues are —1,0,0. We thus conclude that in the 
sense that electric dipoles form an cmtz-ferromagnet (for 
not too oblong samples) 35 the Hamiltonian from Eq. (fT5| 
corresponds to a ferromagnet. (This difference can be 
understood as stemming from the notion that in phonon- 
mediated interactions, it is Zifce-charges that attract, in 
contrast with the Coulomb law.) 

In a mean-field limit, when the interaction J a p(r) from 
Eq. (fT8|) does not depend on the mutual separation r 
between spins, the Hamiltonian (|18[) indeed becomes a 
simple ferromagnet. This limit may be formally imple- 
mented by taking the fc m j cro — > limit in Eq. (|19|) while 
rcnormalizing the integrand so that the total energy of 
the system is finite and scales with the total number of 
the 6-spins. Under these circumstances, the fc-integration 
in Eq. (fT9]) reduces to angular averaging of the K ten- 
sor. Upon the averaging, as outlined in Appendix [B] 
the K tensor becomes a diagonal matrix with eigenval- 
ues that turn out to be dimensionless combinations of 
the microscopic elastic constants: one eigenvalue corre- 
sponding to uniform compression: AT/(A + 2/i), and a five- 
fold degenerate eigenvalue corresponding to pure shear: 
(2/5)(K + 2fi)/(X + 2^i). In other words, in the mean field 
limit, the the dilatational components of local stress on 
one site do not interact with the isochoric components of 
local stress on another site. 

As a result, the mean- field limit of the Hamiltonian 
JH is: 



'compr 

2N 



' N 



T u 



2N 



a=2 



■ N 



,{m) 



(21) 

where s*( m ) = g^/g are unit vectors and denominators 
2Af were introduced for convenience. Note that the cou- 
pling constants J above are proportional to g 2 , according 
to Eq. flTSl) . The numerical value of the </ C ompr/</shear ra- 
tio is equal to the ratio of the two distinct eigenvalues of 
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(a) J > J 

' shear compr 



(b) J . < J 

v ' shear compr 



FIG. 3: Free energy from Eq. (|24[) corresponding to the elas- 
tic constants of (a) silica (Jcompr //shear — 0.88) and (b) sa- 
lol (Jcompr/Jshear — 1-5), corresponding to a high and low 
values of the shear modulus, relative to the bulk modulus. 
Cases (a) and (b) exhibit frozen-in shear and uniform com- 
pression/dilation respectively. Solid circles at bottom plane 
denote the locations of minima. 



the angular-averaged K tensor: 





f 'shear ^/ 






5-component / 
Heisenberg magnet / 


6 












paramagnet 


Ising magnet 





6 / AT 

compr Is 



FIG. 4: The phase diagram corresponding to the mean-field 
Hamiltonian (|21[1 . The shaded region schematically denotes 
the regime where the mapping between liquid dynamics and 
the spin model does not apply. Labels "Ising magnet" and 
"Heisenberg magnet" indicate ordering of the compression and 
shear components in Eq, (12111 in the respective regions of the 
diagram. 



J< 



compr 



K 



_ 5 1 + cr _ 5 3-4c?/cf 
J sh oar ~~ 2 K + 2/x = 2 4-5o- = 2 3 + 2c t 2 /cf 



(22) ing free energy: 



This ratio varies between 0.88 and 1.6 for non-metallic 
glasses surveyed by Novikov and Sokolov. 36,37 Here, a is 
the Poisson ratio, q and Ct are the longitudinal and trans- 
verse speeds of sound. Bigger values of the J com pi-/ «/shear 
ratio correspond to a greater Poisson ratio and smaller 
n/K ratio, i.e. to a lower shear modulus relative to 
the bulk modulus. We thus observe that the mean- 
field Heisenberg model |2TT) has anisotropic couplings, 
the anisotropy directly related to the Poisson ratio of 
the material. Remarkably, the purely isotropic case 
Jcompr/ Jshear = 1 (corresponding to a = 1/5, c 2 /c 2 = 
3/8) falls within the experimental range, the implications 
to be discussed at the end of the article. 

The mean-field Hamiltonian (|2i"j) can be solved in a 
standard fashion by a Hubbard-Stratonovich transforma- 
tion, so that the partition function is given the following 
six-dimensional integral, up to a multiplicative constant: 



Z MF = /n 6 =1 d/i Q exp|-iV 



2/3J Ci 



6, 2 6 N 

a— 2 cy—2 m 



(23) 



where the sum in front of the integral denotes aver- 
aging with respect to spin orientations. (Eq. (|2"3")) is 
a minor variation on the mean field solution of the 
isotropic Heisenberg model for an arbitrary number of 
vector-components, which can be found, for instance, in 
Refj2£) Upon the angular averaging (see the derivation 
of Eq.(fI2]l). the integration can be done by steepest de- 
scent. In the leading order, the partition function is equal 
to e _/3JVG ° , where Go is the minimum value of the follow- 



G 
N 



2(3J C . 



2/3 (/shear 



In 



hi + hi ■ 



(24) 



with respect to the variables h c and h s (c.f. Eq.(2.3) of 
Refj2&). h c = hi gives the effective field on each spin in 
the direction of uniform compression/dilation; it can be 

of either sign. h s = (J2 a =2 ha) 1 ^ 2 1S tne magnitude of the 
total field in the direction of pure shear and can be only 
non- negative, of course. The free energy from Eq. (|24j) is 
graphed in Fig|3] for two distinct values of anisotropy, at 
a temperature below the Curie temperature. 

The Curie temperature itself, according to Eq. ([24|) . is 
determined by the bigger of the two coupling constants: 



7 rri^ 



max(t/ com p r , </ s l lcal ). 

6 



(25) 



We remind the reader that the activated regime corre- 
sponds to the ordered state of the magnet. If J CO mpr > 
•/shear, the ordering is along the compression/dilation 
component of the local displacement, which corresponds 
to the first sum in Eq. ([2"i"T) . see Fig El Note the depths of 
the two minima are equal, implying there is no net vol- 
ume change during the ordering transition. From the spin 
perspective, this state is essentially an ordered Ising fer- 
romagnet. Conversely, if J s hcar > Jcompr, it is the shear 
component that becomes ordered; this component corre- 
sponds to the second sum in Eq. ()21j) . The frozen-in shear 
state is an ordered 5-component Heisenberg ferromagnet, 
from the spin viewpoint. We summarize these notions 
graphically in the phase diagram shown in FigJH 

We may further connect the mean-field parameters 
•/shear and Jcompr to the material constants by enforcing 
the aforementioned notion that the large g asymptotic of 
the function -Fb(ff) be logarithmic in the large g limit, i.e. 
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FIG. 5: The dashed-dotted line shows the mean-field expres- 
sion for the Fo(g) function. The thick solid line shows the 
full free energy from Eq.Q, with a specific choice of the func- 
tion F cx (g)/N = c(g 2 a 3 ), where the coefficient c = 0.05. The 
thin solid line is the common tangent of the two portions on 
the F(g) that correspond to the mechanically stable refer- 
ence state and the aperiodic metastable state. The dashed 
line shows the curve with the largest value of this coefficient 
(c = 5/3), at which the F(g) curve still exhibits a non-concave 
portion. In the inset, we plot the second derivative of the 
mean-field Fo(g). The discontinuity is at g 2 a 3 /ksT = 6, c.f. 
Eq.pTll. 



when the vector length is large. According to Eqs. fH)) 
and the quadratic term in g from Eqs. (ITi|) will can- 
cel out at large g, if 



As a result, 



max(t/ com p r , t/shcar) — 9 Q • 



a _ 6k B TZ ia 



(26) 
(27) 



After comparing this equation with Eq. ([T7j) . we conclude 
that the mean-field Curie temperature of our magnet is 
significantly higher than the temperature at which the 
aperiodic stress pattern sets in, since the actual value of 
g in the metastable phase is significantly greater than 
the critical value from Eq. ([T7|) . see Fig (5] This notion 
confirms that, at least in the mean- field limit, our map- 
ping activated transport onto an ordered spin model is 
internally consistent. 

We further use the relation in Eq. ([26l) to compute the 
mean- field expression for Fo(g) (by numerically minimiz- 
ing the free energy from Eq. ([24l) ). The result is shown 
with the dash-dotted line in FigJSJ To compute the full 
free energy of the solid F{g), from Eq.Q, we must make 
a specific assumption on the penalty F cx (g) for dilating 
the sample. For the sake of argument, we have used a 
quadratic function F ex (g)/Na 3 = eg 2 . A specific value 
c = 0.05 was chosen, so that the resulting barrier be- 
tween the two phases is about fc^T, since it is known 
that the one-particle barrier for surface melting is about 

The mean-field solution supports the asymptotic analy- 
sis from Section [IT] that the fully mechanically stable solid 



is indeed the lowest free energy state; whereas the ape- 
riodic state, though significantly higher in free energy, is 
metastable owing to the surface tension between the two 
phases. 



IV. SUMMARY AND DISCUSSION 

We have mapped the liquid activated-transport regime 
onto the dynamics of a spin model on a fixed lattice. 
The long-lived aperiodic arrangements, which are char- 
acteristic of the activated regime, correspond to aperi- 
odic spin configurations below the Curie point of the 
spin system. The presence of non-zero shear resistance 
in the long-lived structures is reflected in the spins hav- 
ing six components, which correspond to the six inde- 
pendent entries of the local deformation tensor. Each 
spin corresponds to a rigid, weakly interacting molecu- 
lar units (or "bead") of the ROFT theory. The size of 
the bead is set unambiguously as the length scale of the 
inhomogeneity of a frozen-in stress distribution in the ac- 
tivated regime. The length of the vectors is determined 
self-consistently as the value of the order parameter for 
the periodic-to-aperiodic crystal transition. The mutual 
angles between spins change little during the cooperative 
reconfigurations, implying individual bonds are deformed 
nearly harmonically most of the time. The anharmonic- 
ity of the reconfigurations is reflected in a high nucle- 
ation barrier in the free energy profile of the transition, 
the progress coordinate being the number of beads that 
have already switched their positions to their value in the 
new arrangement. The present approach exhibits clear 
parallels with earlier, self-consistent phonon theories of 
aperiodic crystalsji^ which considered the uniform liq- 
uid state as the reference state and have demonstrated 
that the uniform liquid can cross-over to an aperiodic- 
crystal state, accompanied by a discontinuous change of 
an order parameter corresponding to the bead localiza- 
tion length. In contrast, here we construct the aperi- 
odic crystal using a fully mechanically-stable state as the 
starting point, thus allowing for a straightforward treat- 
ment of local, high-frequency shear resistance. We have 
also explicitly demonstrated that aperiodic crystals ex- 
hibit multiply degenerate states whose mutual intercon- 
versions can be mapped onto the dynamics of long-range, 
6-spin model with anisotropic interactions. 

The present microscopic picture has implications for 
the temperature dependence of the activated barriers 
in supercooled melts that we can begin to discuss al- 
ready at the mean-field level. In the allowed part of 
the phase diagram in FigHJ consider a transition from 



the Jcompr < ./shear sector to the J c , 



tor. As already mentioned, this transition would corre- 
spond to ordering of an Ising-like assembly of spins and 
and vice versa for the 5-component Heisenberg spins. To 
avoid confusion, we emphasize that this Ising-like order- 
ing takes place along one component of randomly ori- 
ented, sii-component vectors. Now, according to accu- 
rate calculations^ the heat capacity jump upon order- 
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ing in the classical Heisenberg model tends to increase 
with lowering the dimensionality of the spins and reaches 
its largest value in the Ising model. This implies that the 
r.h.s. part of the phase diagram should generically exhibit 
a higher specific heat. In view of the RFOT-derived con- 
nection between the fragility index and the heat capacity 
jump, i.e. the mentioned m ~ 34.7Ac p relation, we may 
expect that substances with a larger J compI / J s hear ratio 
will generically be more fragile, but under several, rather 
restrictive circumstances, as we discuss next. 

First, we reiterate the present approach is a minimalist 
way to explicitly account for the interactions that give 
rise to (high frequency) shear resistance in supercooled 
liquids. Already in this minimal model, the entropy and 
heat capacity per bead depend on the bead size a, which 
is, formally, the ultraviolet cut-off in the theory. In ad- 
dition to the purely volumetric affect - one vector occu- 
pies a volume a 3 - the bead size will affect the precise, 
self-consistently determined value of the order parameter 
g in the glass phase. Further, the isotropic assumption 
for local elastic response, from Eq. (fTT1) . is an approxi- 
mation. Contributions other than purely elastic terms 
in Eq.© will be generally present as well. One source 
of such contribution is reconfiguration-induced electric 
dipole moment. Earlier estimates 3 ^ suggest that even 
though the local polarization resulting from bead move- 
ments contributes only about a percent to local elastic 
constants, the resulting elemental electric dipoles inter- 
act comparably strongly to the elastic dipoles of the type 
we have considered. This effect should be especially sig- 
nificant in ionic glasses.— Interestingly, because of the 
mentioned antiferromagnetism of the electric dipoles, the 
electric and elastic interactions may be mutually frustrat- 
ing. A careful treatment of polymeric materials, on the 
other hand, should include the effects of chain rigidity 
and other types of local anisotropy that could not be ac- 
counted for by only two elastic constants and the bead 
size a, which are the parameters of our ansatz. Novikov 
and Sokolov have argued electronic contributions must be 
considered in metallic glasses.— 

According to FigJH we should expect that actual sub- 
stances will exhibit a continuous range of local stress dis- 
tributions and heat capacities. Indeed, to distinct points 
on the phase diagram in Fig 01 there correspond very dif- 
ferent values of the heat capacity. Furthermore, the heat 
capacity of an anisotropic Heisenberg model on a peri- 
odic lattice should exhibit a non-monotonic behavior as 
a function of anisotropy, namely a spike at a transition 
that would occur at J CO mpr = Jshear in the mean-field 
limit. Jcompr = ^shcar corresponds to the following value 
of the ratio of the longitudinal and transverse velocities: 
ci/c t ~ 1.63. This is roughly consistent with a very broad 
range of fragilities observed for substances from a rela- 
tively narrow range of the q/c* ratio centered at 1.75 or 
so.— We point out that the correlation between the Pois- 
son ratio and the fragility has been a subject of debate 
for some time i 36 ' 37 i 40 ~ — Our results suggest that a corre- 
lation between the fragility and the Poisson ratio might 
in fact be expected for non-metallic substances that are 



far from the transition region J CO mpr = Jshear- Yet for 
substances near the "critical region" J compl = J s hoar, lit- 
tle correlation is expected. Note that the view of the 
activated regime in a fragile substance as a random Ising 
model below symmetry breaking is consistent with re- 
cent results of Stevenson at el4^ who have mapped the 
localization transition in fiducial liquid structures of a 
Lennard- Jones mixture onto a replica-symmetry break- 
ing transition in a random Ising model. 

Regardless of the detailed value of the heat capacity, 
our results indicate that supercooled liquids exhibit local 
stress distribution ranging from frozen-in compression to 
frozen-in shear. Furthermore, the present findings that 
the activated liquid regime arises self-consistently from a 
mechanically stable reference state merge nicely with the 
self-consistent phonon view of the emergence of the aperi- 
odic crystal state from the uniform UquidJ*^ This notion 
suggests that a unified, quantitative treatment of the liq- 
uid, aperiodic-crystal, and periodic-crystal regimes is in 
sight. We conclude by reiterating that despite system- 
specific variations in direct interactions, accounted for, 
semi-empirically, in the variation of the elastic constants 
in FigUl the underlying mechanism of activated transport 
in non-polymeric liquids is system-independent, consis- 
tent with the conclusions of the RFOT theory. 
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Appendix A: Voigt-Mandel notation 

It is sometimes convenient to present the six indepen- 
dent entries of a symmetric tensor Uij as a six-component 
vector: 

/ "ii \ 

V2u 23 ■ ( Al ) 
V V2ui2 / 



= 1,2,3} i 



so that the convolution h-ijki uu is expressed as bilin- 
ear form Yla,p=i ta^aptpi where A a p is a square 6x6 



matrix: 




(A2) 



2A 3112 
2A 1212 



In the isotropic case, see Eq. ffTT]) , this matrix is rather 
sparse and easily shown to have one eigenvalue 3K and 
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a five-fold degenerate eigenvalue 2/i. Importantly, the 
corresponding orthogonal transformation matrix, i.e the 
matrix consisting of the eigenvectors of the matrix from 
Eq. (|A2[) for an isotropic A, does not depend on the Kj \i 
ratio. (The transformation matrix is easy to compute, 
but is too bulky to give here.) This means, among other 
things, that in the Voigt basis, all isotropic interactions 
are diagonal. 

One may further rescale the components of the 6- 
vectors, for a specific combination of K and /i: 



/ VK(d n + d 22 + d 33 ) \ 
y/Jt(du - d 22 ) 
f (dii +d 22 -2d 33 ) 

2^/7^23 
2-^/7^31 

2y/Jid 1 2 



(A3) 



so as to make the matrix from Eq. (IA2[) not only diagonal, 



but also unit: dijAij k id k i = 



E«=i 9l = 3 1 (and 



E ra a 3 d nj (r m ) A nj i p (-iki)u p t 



-ikr n 



where an appropriate ultraviolet cut-off at fcmicro, as in 
Eq.([ll]), is understood. (Note A ijk i = A jiku A ijk i = 
Aijik-) For the specific, simple form of the elastic tensor A 
from Eq. (fTTj) , it is straightforward to compute the Gaus- 
sian integral above explicitly, using the obvious property 
of the Fourier component of the real-valued displace- 
ments u: u(—k) = u*(k), and Hubbard-Stratonovich 
formulas ^ exp[-(a; 2 + y 2 ) + a(x ~ iy) + b*(x + 
iy)] = e ab " and J^ ocj (U m dx m /T: 1 ^)e- XkAklXl+bkXk = 

(DetA)- 1 / 2 e bfcj4 ^ lb,/4 , where {A kl \ is a positively de- 
fined symmetric matrix. The inverse of the matrix An = 



Aijimkjkm is easy to find, using the special property of 
the matrix Pj m = kjk m /k 2 = k jk m that P 2 = P. Using 
this notion, one may show that 



(Aijlm kj k n 



(B3) 



where Aijki is from Eq. tfTTj) . Finally note that the 
fc-integration in the exponent in Eq. (|B2j) counts every 
phononic mode twice, because of the mentioned property 
u(-k) = u*{k). (d 3 u(k) = {d 3 Re[u(k)}d 3 lm[u{k)}}^ 2 .) 
We use the listed notions to perform the functional inte- 
gration in Eq. (|B2|) and thus arrive at Eq. fp~8|) of the main 
text. 

It may be instructive to compare the straightforward, 
if not somewhat tedious calculation above with a simpler 
case of interaction between 3-component vectors medi- 
ated by a single-polarization elastic interaction: 



e = y"di/[ J ft:(vv) 2 /2 + Xl(f (m)v?/ ') j3 ( r " 



.)]» I 54 ) 



similarly for the elastic components), see Eqs.([§] [TU|) . 



Appendix B: Spin-spin interaction 

To evaluate the coupling between the non-linear oscil- 
lators, we integrate out only the phonon degrees of free- 
dom u in Eq. (J4j , i.e. we need to perform the following 
functional integral: 

'H r d 3 u(r)]e~ 13 ^ d3r ^ jAzjk ' Uk ^ 2 ^ + ^ m a 3 dij(r m )A ijk iu M {r m )] 

(Bl) 

where the m-summation is over the locations r m of 
the harmonic oscillators. (Note the dAu cross-term is 
symmetrized because Ay/.; = A k uj.) In terms of the 
Fourier components of the local displacements u: u{k) = 
J d 3 r u(r)e ikr , this integral becomes: 

ju k d 3 u(k)}exp^-l3 J ^3[A nj i m u n (k)ui(-k)kjk m /2 



where ip = tp(r) is a coordinate dependent, scalar field. 
In this case, the phonon-mediated coupling between the 
spins becomes: 



n 



^ ] a Jij{v rnn )g i 



'-) in) 



(B5) 



This coupling looks particularly simple in the Fourier do- 
main: Jij(r mn ) = J T2~j3 cos (fcr ) ki kj ; at distances well- 
exceeding the inverse ultraviolet cut-off, it is identical to 
, the usual electric dipole-dipole interaction, albeit with 
the opposite sign: Jij(r mn ) = (% - 3mrij) /r^, where 
r mn = (r m - r„) and n = r rnn /r mn . It is easy to show, 
by a direct calculation, that the eigenvalues of the matrix 
kikj are 1, 0, 0. 

To perform the angular averaging of the tensor K in 
Eq. (|19[) . one needs to compute {kikj) and (kikjkik m ), 
where the angular brackets (...) denote angular averag- 
ing. To compute these averages, we note that by sym- 



]}, (B2) 



metry, (kikj 



oc 5ij and that kf 



1, implying (kikj 



(l/3)6ij. By the same token, since (kikjklk m ) OC (o~ijO~Lrn~\~~ 
SilSjm + S im 5ji) and k 2 k 2 = 1, we obtain (kikjhk m ) = 
(1 / '15) (Si j Sim + SuSjm + SimSji). Using these expressions 
and Eq.lJTTJl, we obtain (Kijkl) = XSijSu + (J,'(6ik5ji + 
SuSj k ), where A' = (A 2 + 16AaV15-V/15)/(A + 2^) and 
H 1 = 2/i(3A+8/x)/15(A+2/i). Therefore the (K ijH ) tensor 
has the same structure as the A tensor itself and will be 
diagonal in the Voigt representation, after a basis change 
as in Eq. (|A3|) . Of the corresponding diagonal entries, one 
is K/(\+2(i) and the rest five are (2/5)(# +2ju)/(A+2^), 
as mentioned in the main text. We note that the sum of 
the diagonal entries is equal to 3, as expected from the 
discussion in the beginning of Section [ 
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